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Abstract — The achievable data rates of current fiber-optic 
wavelength-division-multiplexing (WDM) systems are limited by 
nonlinear interactions between different subchannels. Recently, it 
was thus proposed to replace the conventional Fourier transform 
in WDM systems with an appropriately defined nonlinear Fourier 
transform (NFT). The computational complexity of NFTs is a 
topic of current research. In this paper, a fast inverse NFT 
algorithm for the important special case of multi-solitonic signals 
is presented. The algorithm requires only 0(3} log 2 37) floating 
point operations to compute D samples of a multi-soliton. To the 
best of our knowledge, this is the first algorithm for this problem 
with log 2 -linear complexity. The paper also includes a many 
samples analysis of the generated nonlinear Fourier spectra. 

Index Terms — Nonlinear Fourier Transform, Fast Algorithms, 
Optical fiber, Nonlinear Schrodinger Equation, Solitons 

I. Introduction 

The current generation of fiber-optic communication sys¬ 
tems based on wavelength-division-multiplexing (WDM) is 
operating close to capacity m. In order to avoid a looming 
capacity crunch, new approaches beyond WDM need to be 
explored. The nonlinear effects in optical fiber lead to interfer¬ 
ence between different subchannels, and have been identified 
as one of the major capacity-limiting factors. In order to avoid 
interference between subchannels, it was recently proposed 
to modulate information using a nonlinear Fourier transform 
(NFT) instead of the traditional Fourier transform used in 
WDM [|2|, [ED, 0. Analogously to how the spatial evolution 
of a signal in a linear channel such as copper cable becomes 
trivial in the traditional Fourier domain, the evolution of a 
signal in optical fiber becomes trivial in a properly defined 
nonlinear Fourier domain. The most important aspect here 
is that different regions in the nonlinear Fourier domain do 
not interact as the signal travels through the fiber. Thus, 
by defining subchannels in the nonlinear Fourier domain, 
interference among them can be avoided. There are many 
interesting open questions regarding information transmission 
in the nonlinear Fourier domain, with one of them being 
computational complexity. While the celebrated fast Fourier 
transform has been the workhorse of practical implementations 
of WDM, the first fast NFTs have been proposed only recently 
0,0- The inverse NFT however, in contrast to the traditional 
case, differs from the forward transform. No fast inverse NFT 
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has been established in the literature so far to the best of 
our knowledge, although inverse NFTs are used in several 
practical problems f7), J8j]. The goal of this paper is to, 
at least partially, close this gap for the important case of 
multi-solitonic signals. The algorithm developed in this paper 
can generate multi-solitons with arbitrary discrete spectra 
(but not norming constants). Such signals have recently been 
considered for information transmission in 0 and 0. 

The paper is structured as follows. In Section [II] the 
continuous-time NFT is introduced and discretized. The pro¬ 
posed algorithm will operate in two stages, which are then 
developed in the Sections [HI] and [IV] Two numerical examples 
follow in Section |V] Section |VT| finally concludes the paper. 

II. Nonlinear Fourier Transform 

In this section, the nonlinear Fourier transform is first 
established for continuous-time signals. The NFT is then 
discretized and, finally, the inverse discrete NFT is discussed. 


A. Continents-Time Transform 

The nonlinear Fourier transform [w.r.t. the signal model 
in Eq. ([8]) below] is due to Zakharov and Shabat fTOl . Fet 
{q(t)}tfs k <= C denote a vanishing signal. Its NFT is found 
through the analysis of the associated scattering problem 
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With </>i and </>2 denoting the components of <fi, define 


a(A) := lim e lAt 0i(7, A), /3(A) := lim e lAt 02(7, A), (3) 

t —>00 t —>00 

where A e C. The functions a and /3 are usually not well- 
defined on the whole complex plane. With some assumptions 
on the signal q(t), however, the following can be said. 

Lemma 1 (liTTI). If there exist positive constants C and k 
such that |q(7)| < C e~ 2K ^ for all real t, then: 

1) a(A) is analytic in the closed upper half-plane 5(A) ^ 0 

2) a(A) has finitely many roots Ai,..., A j in 5(A) ^ 0 

3) /3(A) is analytic around the real line R. and well-defined 
at the roots Ai,..., A j of a( A) in 5(A) > 0 





This lemma implies that the reflection coefficient 

«( w ) : = weR, (4) 

a(ui J 

is well-defined except at finitely many poles. The function (3/a 
is also well-defined at the (complex) roots A j. The residual of 
/ 3/a at these points is, assuming that all A ; are simple roots, 
given by 

With these definitions in mind, the nonlinear Fourier trans¬ 
form (or scattering data ) of q(t) is defined as ifTTl Apdx. 5] 

(**>&)/=!■ (5) 

Please note that the nonlinear Fourier spectrum has a physical 
interpretation DU. The reflection coefficient reduces to the 
linear Fourier transform for signals with small amplitudes. It 
describes the so-called radiation components of the signal. 
The roots A*, indicate solitonic components in the signal, 
which are further specified through the norming constants 
< 7 fc. (Solitons are basically “travelling humps” that have no 
analogue in linear evolution equations.) A multi-soliton is a 
signal that contains several solitonic components. A signal is 
reflectionless if the reflection coefficient is zero everywhere. 

The nonlinear Fourier transform has several properties sim¬ 
ilar to the conventional Fourier transform |2] Sec. IV.D], In 
particular, the effects of shifting or dilating the signal on the 
nonlinear Fourier transform turn out to be simple: 

q(t-t 0 ) <-> {e“ 2l “‘°g(w)j, qf), (6) 

q(t/v) <->■ {\v\q(vu)}, iv^j, \vWj)- ( 7 ) 


The most important property is however the following. Con¬ 
sider a space-time signal {Q(x,t)} x ^ 0 ,teR c C that is gov¬ 
erned by the focusing nonlinear Schrodinger equation (NSE) 



d 2 Q 
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( 8 ) 


Denote the functions a(z) and (3{z) that correspond to the 
signal q(t) = Q(x,t) by a(x, A) and (3(x, A). Then, |j2] V] 


a(x, A) = a(0, A), (3(x, A) = e 4l A x (3(0, A). 


The focusing NSE describes the evolution of a signal in 
lossless, noise-free optical fiber (3- The evolution of a signal 
thus indeed becomes trivial in the nonlinear Fourier domain. 


B. Discrete-Time Transform 

In order to compute a numerical approximation of the 
nonlinear Fourier spectrum 0, the signal q(t) is considered 
only on a (sufficiently large) finite time interval. Since the 
effects of time-shifting and dilation are simple [see 
we can choose the interval [—1,0] in order to simplify the 
following exposition without loss of generality. In this interval, 
D equidistant scaled samples 

Q[n\ := eq(-l + ne- 0 , e := n e {1,..., D}, 


are taken in order to represent to signal. It will be convenient 
to transform the spectral parameter A as follows: 


z = z( A) := e 
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(The real part of A has been constrained so that the inverse 
transform z X(z) is well-defined.) Under the assumption 
that q[t) is piecewise constant and that e is infinitesimal, the 
scattering problem 0-0 reduces to ]8] Sec. II] 
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Here, the initial condition ( [TT| was obtained by evaluating 
0 at t = — 1. Similarly, after considering 0 at t = 0, the 
functions a(A) and /3(A) are approximated by 


a(z) '■= <t>i[D,z\, b(z) := <j>i[D,z]. (12) 


Note that the complex half-strip 5ft(A) 6 [— ^], 3(A) > 

0 with respect to the original spectral parameter A corresponds 
to the exterior of the unit circle 1*1 > 1. The real interval 
[-■%-, f-], in particular, is mapped to the unit circle \z\ = 1. 
The reflection coefficient q(ui) is therefore approximated as 


Q(u) 


b(z(tu)) r 7T 7T 

a(z(uj)) ’ l 2e ’ 2e 


(13) 


The roots zi, ... ,Zk of a(z) outside the unit circle \z\ > 1 
are taken as approximations of the roots Ai,..., Aj of a(A) 
subject to the coordinate transform 0. Note that there may be 
spurious roots Zk that do no correspond to any of the A j. These 
spurious roots are numerical artifacts. They usually occur in 
physically unplausible regions (e.g., very close to or very far 
away from the unit circle) so that they can be easily filtered 
out. The norming constants q 3 are approximated through 


Qk • b(z k ) j 


da dz . , 
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(14) 


The discrete nonlinear Fourier transform is then defined by 


{Q( w )}o;e[— { z kiQk}k= 


(15) 


C. Inverse Discrete-Time Transform 

A close look at the finite difference equation ( fi~0| ) reveals 
that a(z) and b(z ) defined in (12 1 are neccessarily polynomials 
of z _1 with a degree of at most D |8] (19)]: 


(*) = b{z) = YX=obi 


D—li -i 


(16) 


Please note that the additional factors z ,)/2 given in |8] (19)] 


cancel with the factor in (11 1 in our case. The following lemma 
specifies the range of the discrete nonlinear Fourier transform. 

Lemma 2 (181). Fix two arbirary polynomials as in ( |76| ). Then, 
there exists samples Q[ 1],..., Q[i9] such that (12) is satisfied 















if and only a,Q ^ 0 and |a(£)| 2 + |6(£)| 2 = 1 far all |£| = 1. 
In that case, ag = 1 / n n =i V 1 + iqwi 2 > ° and 
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Proof: The first part in shown in |^1 Sec. IV+Apdx.]. The 
other two formulas are given in © Eqs. (14a)+(15)+(37)]. The 
last formula deviates from © because here Q[n] = —p n - ■ 

There are various ways to approach the inverse nonlinear 
Fourier transform. See, e.g., @. The fast 0(77 log 2 D) algo¬ 
rithm that will be proposed in the following two sections uses 
a two-step procedure similar to © and (121 : 

1) Synthesize a(z) and b(z) such that the conditions in 
Lemma[2]are met and the roots of a(z ) are z k '■= z(X k ). 

2) Recover the samples Q[n] from (10 1 —(11 1 . 


The complexities of comparable algorithms are higher: either 
0(D 2 ), as in ED and ®; 0(DK 2 ) as in © and JU III.B-3]; 
or even 0(D 3 ) as in © III.B-1], Here K is the number of 
solitons. (Note that K ~ D in communication problems.) 


III. Step 1: Fast Generation of a(z) and b(z) 

In this section, the first step of our proposed algorithm is 
described. The basic approach used here is due to our previous 
work (321. However, the construction given here avoids a 
pathology from which the construction in 1321 suffers, as a 
new many-samples analysis presented below shows. Please 
note that no complexity analysis of this step will be given 
here because the analysis given in lfl2l VI] still applies. 


A. Basic Construction 

Consider any polynomial f){z) = Yn=o that fulfills 


is not too small. We shall assume this in the following, and 
set 


a{z) := Pi z 1 = e 1 v a ide ai(-z) + 0(e 2 ), 


where ip e [—7T, 7 r] is such that a(z) fulfills the condition 
lim^oo a{z) ^ 0 in Femma[2] In light of (20 1 , one finds that 
a{z) and h(z) also fulfill the other condition in Femma[2](up 
to a rapidly vanishing error). Thus, the problem of recovering 
<3[1],..., Q[D\ through < 10 >—< 11 > is well-defined. 


B. Choice Of The Filter i/f z) 

In our previous paper lfl~2l . the filter rp(z) = z ( D ^ and a 
very small 5 were used. However, an analysis similar to the one 
below shows that this choice results in a reflection coefficient 
that does not vanish in the many samples regime for w —* ±oo. 
Parseval’s identity for the NFT J2] p. 4320] however implies 
that any relection coefficent of a signal with finite energy must 
vanish. Although this can be fixed by letting S depend on e, one 
then obtains the problem that the norming constants converge 
towards zero as D —* oo. Again, this cannot happen for any 
real signal. In order to avoid these two problems, one may 
fix S independently of e and use a proper low-pass filter for 
4>(z) to ensure that the reflection coefficient vanishes at higher 
frequencies. In our numerical experiments, a simple truncated 
ideal low-pass was used 131 P- 98]: 
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( 21 ) 


where ui c = 10 is the desired cut-off frequency for the 
reflection coefficient ([4]!. Please note that the filter 'ib(z) is 
independent of q(t) and thus needs to be designed only once. 


WOI e [0,1] V|z| = 1. 


Through spectral factorization 03, two polynomials u(z) = 
UiZ~ l and b(z) = biZ~ l can be found that have 

no roots outside the unit circle and satisfiy 


\ U (o\ 2 = i-6 2 \m\ 2 ,m\ 2 = 5 2 \m\ 2 vki 

Here, S e (0,1) is a free parameter. We choose 
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This ideal aideai(z) has the right roots in M > 1 and satisfies 


hdeal(£)| 2 + IKOI 2 =H0\ 2 I! \Z\ l -== + IKOI 2 

fc=l 

(use ICI = 1) =l-5 2 \m\ 2 + S 2 m)\ 2 = h (20) 


for all |£| = 1, but it is not a polynomial in z~ 3 , in general. 

However, since ai ( | ca ] (z) is rational with all poles inside the 
unit circle, it has an expansion a ic i ea i(a:) = Yi=o a ideai,i z ~ l 
in \z\ ^ 1 with coefficients that vanish exponentially vast. 
Thus, Oideai(-) can be easily approximated by a polynomial 
p(z ) = Yb!i=o Pi z ~ l U P to an y desired precision given that D 


C. Many Samples Analysis 

In this subsection, the continous-time nonlinear Fourier 
spectrum |5]) that is generated through a(z) and b(z) in¬ 
troduced above is analyzed in the many samples regime 
D = 1/e —> oo. The function a(z) is constructed such that 
its roots in |z| > 1 match exactly the prespecified Ai,..., A k 
under the transform A z(X). The behavior of the reflection 
coefficient can be analyzed using ( fi~3] i, ( p~8] > and ( [20] ): 

<5 2 |vE-H| 2 
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( 22 ) 


denotes the asymptotic amplitude of (21 1 given by 
= lim |'0(z(w))| = 


Si(w + ui c ) — Si(w — uj c 


D—> oo 


where Si(s) := m j d -dd is the sine integral function lfl4l p. 
98]. It now follows from |'P(w)| —» ±\ that q(fa) —» 0 for 
u> —* ±oo. That is, the pathology of the construction in lfl2l 
discussed in the previous subsection has been fixed with the 
new design. The norming constants (14 1 are more difficult to 
analyze. In the appendix, they are shown to satisfy 


q k = lim Q k = 


B( A fe ) 
U(X k ) 


(Afe — A, 
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where B( A) := lim c ^ 0 b(z(X)) an d (7(A) := lim e _> 0 u(z(A)). 
This formula shows that the norming constants have finite 
non-zero limits in the many samples regime unless one of 
the prespecified A*, coincides with a root of either B( Afc) 
and U( Afc). We have not yet performend an exact analysis 
of B( A) and 17(A). But since the functions B( A) and 17(A) 
(and therefore their roots) are independent of the Ai,..., A k, 
it seems unlikely that the pathology of vanishing norming 
constants discussed above occurs with the new construction. 


IV. Step 2: Fast Inverse Scattering 
In this section, a(z) and b(z) are now considered known 
and the samples Q[l], ..., Q(D] have to be recovered. A 
fast inverse scattering algorithm that recovers Q[ 1],..., Q[D] 
using only 0(D\og 2 D) floating point operations (flops ) for 
systems of the form ( [TO) ) has been proposed by McClary on 
(2)]. Unfortunately, his algorithm does not work in our setup. 
In || 151 , the functions a(z) and b(z) are interpreted as z- 
transforms of a down-going wave, which is sent into the earth, 
and a resulting upcoming wave, which is measured at the same 
location. The initial condition ( fTT) corresponds to the situation 
where a unit pulse is sent into the earth and no upcoming wave 
is generated, leading to a trivial medium that does not reflect 
(i.e., Q[ 1] = • • • = Q[D\ = 0). 

In the following, we will describe how McClary’s approach 
can be modified such that it works in our setup. By inverting 
( 101 and taking the definition ( 121 of a(z) and b(z) into 
account, one arrives at the following scattering problem: 
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Assume that D is even and split 


i(z) = a L (z) + z 




(z), 


(24) 

6[n\ := V 1 + IQMI 2 - 
\z) := YIJo 1 a i z~\ 


Split b(z) = b L (z) + z 2 b (z) in the same way, and define 


T„(*) : = 


V 1 + IQMI 2 
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zQ[n\ 


-Q[n\ 
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(25) 


i{z) :=T„ +1 ( 2 )T n+2 ( 2 ) x • ■ ■ x T m (z). 


With this notation, one has that 

= t d 


OD-r, 
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(26) 


This equation leads to the important conclusion that the highest 
coefficient of ao-m and &D_ m is independent of ajj and bjj 
as long as to < y (because z _ TTo_, m will always shift the 
coefficients of ajj and bjj up in that case.) Since Q[D — to] 
can be recovered from the highest coefficients through ( | 1 1\ , 
this shows that the iteration ( |24| ) with initial condition 

an ( 2 ) = a u (z), bn ( 2 ) = b u (z) 

can be used to recover Q[D],.. ., Q\® + 1]. This insight 
ensures that the divide-and-conquer strategy given in Algo¬ 
rithm [I] indeed recovers Q\ 1],..., Q[D], The computational 


Algorithm 1 Recursive inverse scattering (requires D = 2") 


Input: a(z) = 0 , 2 “* and b(z ) = 

Output: Q[l],... ,Q[D] and Td^i{z) 
if D = 1: 


D — l 


• return Q[l] = ao/bo and Ti( 2 ) in (25 i 
else: 

• call self with inputs a u (z ) and b u (z) to recover 
Q[D], ..., Q[y + 1] and T D ^ D / 2 {z) 

• recover a D / 2 (z) and b D / 2 {z) via ( 261 

. call self with inputs ar>/ 2 (z ) ando 0 / 2 ( 2 ) to recover 
Q[D/2],..., Q[l] and T D/2 ^ 

. return Q[l],... ,Q[D] and T D ^ D / 2 T D / 2 ^\ 


bottleneck in Algorithm [T] are the polynomial products. By 
applying the same complexity analysis that was used for El 
Alg. 1 ] to the recursion tree generated by Algorithm |T| one 
finds that its overall complexity is 0(D log 2 D) flops if the fast 
Fourier transform is used to multiply polynomials efficiently. 

V. Numerical Examples 

In the first example, a single soliton with a root at Ai = 20 i 
is considered. For small S > 0, the generated signal will be 
almost reflectionless. For several number of samples D and 
S’ s, the two steps of our algorithm as described in Sections 
uni and [TV] have been carried out. Then, the forward NFT was 
used to compute the norming constants of the generated signal. 
Reflectionless signals with one soliton have a well-known 
closed-form solution (e.g., fl2l II.D]). We used this analytical 
solution to generate samples Qcxact [ n \ that represent the exact 
one-soliton with the prespecified Ai and the numerically 
determined norming constant. The error between the generated 
signal and the analytical solution is depicted in Figure [T] (left). 
It eventually saturates because the generated signal is not 
exactly reflectionless, but the error floors decrease as 5 is 
decreased. The normalized runtimes are presented in Figure [I] 
(right). The runtime per sample grow only very slowly with D, 
confirming the overall 0(1) log 2 D) runtime. (The runtimes do 
not include the generation of u(z) and B(z) because they are 
independent of the signal and have to be computed only once.) 

The second example considers a problem from 0. A 
multi-soliton with roots at A^ = 0.25 i fc, fee {1, ..., 4}, is 
generated. The roots were scaled by a factor 77 = 50 such 
that the signal would fit in the interval [—1,0] (see (|7]i). The 
generated signal was rescaled to the original time-scale and 
analyzed using the forward NFT. The results are depicted in 
Figure [2] The forward NFT found the roots Afc of the generated 
signal exactly in the right places. The reflection coefficient 
looks as predicted in ( [22| . It is interesting to compared Figure 
[2] (left) with [9] Fig. 2b], It was expected that the signals 
are different due to different norming constants. However, our 
signal has approximately the same duration, but with only half 
as many extrema. The highest peak is also about only half as 
high as in J9] Fig. 2b], We speculate that the roots in (9] 
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Figure 1. First example. Left: Squared error £n=i IQM — Qexact[rr]| 2 /2 |Qexact[n]| 2 - Right: Runtime per sample. 




Figure 2. Second example. Left: Generated time-domain signal (D = 512, <5 = 0.01). Right: Roots of a(z(X)) and reflection coefficient. 


Fig. 2b] may be double because these quantitative differences 
vanish if the orders of the roots are doubled in our algorithm. 

VI. Conclusions 

The (to the best of our knowledge) first fast inverse non¬ 
linear Fourier transform for the generation of multi-solitons 
has been presented. The generated nonlinear Fourier spectra 
have been analyzed through a many samples analysis. The 
0(D log 2 D) runtime and accuracy of the algorithm have also 
been confirmend experimentally. The current algorithm offers 
only rough control over the reflection coefficient and norming 
constants, which should be addressed in future research. 
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Appendix: Derivation Of Equation ( |23| ) 

The transformed coordinate z = e _2lAc will be close to 
one in the many samples regime e = 1/D —* 0. The usual 
first order approximation of the exponential around zero, 

°° J 

eS = ^]Tp = l + s + 0(|s|"), 
i =o *■ 

shows that 

z = 1 — 2i A e + 0(e 2 ) 

and 

z = e 2i * e = 1 + 2 i A e + 0(e 2 ). 


From these two formulas, one finds that 

z — Zi (1 — 2 i Ae + 0{e 2 )) — (1 — 2 i \e + 0(e 2 )) 
1 -zzi = 1 - (1 - 2 i Ae + 0(e 2 ))(l + 2 i A iC + C>(e 2 )) 
A — Xi + 0(e 2 ) 

~ X-K + 0(e 2 ) 
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1 - z k z k 1 - (1 - 2i A fc e + 0(e 2 ))( 1 + 2i A fc e + 0(e 2 )) 
1 1 


2 ie Afc — Afe + 0{e 2 ) 

The derivative of the Blaschke product 

B{z) := [] 


Z-Zj 
1 - ZZi 


that occurs in the definition < [l9| i of a;d e ai(z) is given by 

dB . N _ 1 — \z k \ 2 -p-r z — Zi 


and satisfies 


dB 

dz 


(Zk) 


1 

i-k-l 2 


n 

i^k 


Zk ~ Zj 
1 - Z k Zi ' 


The approximations given above therefore imply that 
dB, , _ 1 1 -p-r Afe — A; + 0(e 2 ) 

dz [Zk) = 2ie A fc - A fc + 0(e 2 ) H ~ X k - X, + 0(e 2 ) ’ 

and thus 

^(z k )= da f^(z k ) + 0(e 2 ) 
dz dz 

= ^(z k ) B{z k ) +u(z k )~(z k ) + 0(e 2 ) 
dz ^ —v—' dz 

=o 

_u(z k ) 1 -p-r X k — X i+0(e 2 ) 

- 2 i e A fc - A fe + 0(e 2 ) H X k - A, : + 0(e 2 ) 

+ 0(e 2 ). 
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Figure 3. Norming constants for the second numerical experiment 


n 

iy£k 


Afc — A i + 0(e 2 ) 
Afc - Ai + 0(e 2 ) 


b(z k ) / (u(z k ) 


z k 


+ 0(e 


1 


1 Afc - Afc + 0(e 2 ) 


n- 

i^k 


Afc — Xi + 0(e 2 ) 
Afc - A* + 0(e 2 ) 


B( Afc) 
U( Afc) 


(Afc 


^fc) n 


+ 0(e 2 )) 

X k Ai 

Afc Ai 


where U{ A) and B( A) have been defined in Section III-C 


We have not yet analyzed the limits U{ A) and B( A), as was 
already mentioned in Section [TlI-C| However, we compared the 
prediction made by the semi-asymptotic formula 




u(z k ) 


i^k 


(27) 


with the exact norming constants that were determined nu¬ 
merically for the second numerical experiment in order to 
investigate how fast the norming constants converge towards 
their many samples limit. The results are shown in Figure 
0 The values obtained from <[27]) match the true norming 
constants quite well, but not as well as the roots match their 
desired values in Figure [2] (right). This is consistent with the 
results reported in ||3] p. 4339], were it was found that the 
norming constants (or spectral amplitudes ) converge slowly. 


Equation ( |23j ) now follows when the previous formula is 
plugged into (|14|>: 


Qk 


b(z k ) 


u(z k )_1_ 

2 i e Afc — Afc + 0(e 2 ) 


2 i ez k 






































